calcSummaryOutputValues_ <- function(x, t, n, useXRangeOffset=F) {
  
  if (!(useXRangeOffset)) {
    l=min(x)
    u=max(x)
  } else {
    # the offset here is to handle division by zero
    l=min(x)+0.00001  # lower end of x range where we assume linear contextual model
    u=max(x)-0.00001  # upper end of x range where we assume linear contextual model
  }
  
  p=length(x) # number of precincts
  
  nx1=sum(n*x)
  ###########################
  # DD bounds:
  bl <- numeric(p)
  bu <- numeric(p)
  for (i in 1:p) {
    lowerBound <- max( 0, (t[i] + x[i] - 1) / x[i] )
    upperBound <- min(1, t[i] / x[i] )
    if (!is.finite(lowerBound)) {
      lowerBound <- 0.0
    }
    if (!is.finite(upperBound)) {
      upperBound <- 1.0
    }
    bl[i] <- lowerBound
    bu[i] <- upperBound
  }
  ###########################
  # mean of midpoint of deterministic bound (for reference); note that this is not scaled/weighted by the population
  ddMidPointMean <- mean((bu-bl)/2 + bl)
  ###########################
  
  r=sum(n*x*(1-x))/sum(n*x)
  lb=0 # this is our choice of lambda (always 0 in the current version)
  h0=lb*sum(n*x*t)/sum(n*x)
  s1= sqrt(sum((n*x*((1+lb)/2-lb*x))^2)/(sum(n*x))^2)
  cl=0.90
  #this is conf level
  ex=1
  gl01=0
  gl1=c(-1/l,0,0)
  gl02=-1/(1-l)
  gl2=c(1/(1-l),1/(1-l),1/(1-l)-1)
  gu01=1/l
  gu1=c(-1/l,0,0)
  gu02=0
  gu2=c(1/(1-l),1/(1-l),1/(1-l)-1)
  ###
  gl03=0
  gl3=c(-1/u,0,0)
  gl04=-1/(1-u)
  gl4=c(1/(1-u),1/(1-u),1/(1-u)-1)
  gu03=1/u
  gu3=c(-1/u,0,0)
  gu04=0
  gu4=c(1/(1-u),1/(1-u),1/(1-u)-1)
  h=c(1-lb, sum(n*x*(1-lb*x))/sum(n*x),sum(n*x^2*(1-lb*x))/sum(n*x))
  ######
  
  #in practice w0 should be estimated by regressing t<-w0+c1x+d1x^2
  #d1=b1-w1 #in practice d1 should be estimated by regressing t<-w0+c1x+d1x^2
  #c1=b0-w0+w1  #in practice c1 should be estimated by regressing t<-w0+c1x+d1x^2
  reg=lm(formula = t ~ poly(x, 2,raw=TRUE)) # do quadratic regression here
  w0=coef(reg)[1]
  c1=coef(reg)[2]
  d1=coef(reg)[3]
  
  th=c(w0,c1,d1)
  v=sandwich::vcovHC(reg,"HC1")
  #this gives sandwich variance of STATA see
  #https://stats.stackexchange.com/questions/117052/replicating-statas-robust-option-in-r
  #information re robust sandwich variance formula:
  #https://www.stata.com/manuals/p_robust.pdf
  #
  sl1=s1+sqrt( t(h-r*gu1)%*%v%*% (h-r*gu1))
  sl2=s1+sqrt( t(h-r*gu2)%*%v%*% (h-r*gu2))
  sl3=s1+sqrt( t(h-r*gu3)%*%v%*% (h-r*gu3))
  sl4=s1+sqrt( t(h-r*gu4)%*%v%*% (h-r*gu4))
  sl=c(sl1,sl2,sl3,sl4)
  su1=s1+sqrt( t(h-r*gl1)%*%v%*% (h-r*gl1))
  su2=s1+sqrt( t(h-r*gl2)%*%v%*% (h-r*gl2))
  su3=s1+sqrt( t(h-r*gl3)%*%v%*% (h-r*gl3))
  su4=s1+sqrt( t(h-r*gl4)%*%v%*% (h-r*gl4))
  su=c(su1,su2,su3,su4)
  bl1=h0-r*gu01+t(h-r*gu1)%*%th
  bl2=h0-r*gu02+t(h-r*gu2)%*%th
  bl3=h0-r*gu03+t(h-r*gu3)%*%th
  bl4=h0-r*gu04+t(h-r*gu4)%*%th
  bbl=c(bl1,bl2,bl3,bl4)
  bu1=h0-r*gl01+t(h-r*gl1)%*%th
  bu2=h0-r*gl02+t(h-r*gl2)%*%th
  bu3=h0-r*gl03+t(h-r*gl3)%*%th
  bu4=h0-r*gl04+t(h-r*gl4)%*%th
  bbu=c(bu1,bu2,bu3,bu4)
  
  wuc=c(gu01+t(gu1)%*%th,gu02+t(gu2)%*%th,gu03+t(gu3)%*%th,gu04+t(gu4)%*%th)
  wu=min(gu01+t(gu1)%*%th,gu02+t(gu2)%*%th,gu03+t(gu3)%*%th,gu04+t(gu4)%*%th)
  #upperbound of w1
  wl=max(gl01+t(gl1)%*%th,gl02+t(gl2)%*%th,gl03+t(gl3)%*%th,gl04+t(gl4)%*%th)
  #lowerbound of w1
  
  cil=max(bbl)-ex*sl[which.max(bbl)]
  cir=min(bbu)+ex*su[which.min(bbu)]
  
  #these are the conservative ci at ex=1.
  ###########################
  
  hbdu0=min(bbu) #hat district level upperbound proposed (wtd)
  hbdl0=max(bbl) #hat district level lowererbound proposed (wtd)
  
  bdu=sum(n*x*bu)/sum(n*x) # district level upperbound by duncan-davis (wtd)
  bdl=sum(n*x*bl)/sum(n*x) # district level upperbound by duncan-davis (wtd)
  
  outputList <- list()
  outputList[["nx1"]] <- nx1 # total elements (people/households/etc.) of variable X across all geographic units
  outputList[["hbdl0"]] <- hbdl0 # CI_0 lower
  outputList[["hbdu0"]] <- hbdu0 # CI_0 upper
  outputList[["cil"]] <- cil # CI_1 lower
  outputList[["cir"]] <- cir # CI_1 upper
  outputList[["bdl"]] <- bdl # DD lower
  outputList[["bdu"]] <- bdu ## DD upper
  # for reference/analysis:
  outputList[["w0"]] <- w0
  outputList[["c1"]] <- c1
  outputList[["d1"]] <- d1
  outputList[["p"]] <- p
  outputList[["l"]] <- l
  outputList[["u"]] <- u
  outputList[['w1_l']] <- wl
  outputList[['w1_u']] <- wu
  outputList[['b0_u']] <- c1 + w0 - wu
  outputList[['b0_l']] <- c1 + w0 - wl
  outputList[['b1_u']] <- d1 + wu
  outputList[['b1_l']] <- d1 + wl
  outputList[["ddMidPointMean"]] <- ddMidPointMean
  
  return(outputList)
  
  
}



ei_bounds = function(x, t0, t1, n, k = 100) {
  
  #get diff-in-diff estimate
  m_did = lm(t1 - t0 ~ x, weights = n)
  att = m_did$coefficients[2]
  #Get king et al bounds
  bounds_t0 = calcSummaryOutputValues_(x, t0, n)
  bounds_t1 =  calcSummaryOutputValues_(x, t1, n)
  
  #Create grids
  l = min(x)
  u = max(x)
  
  seq_w1_t0 = seq(bounds_t0$w1_l, bounds_t0$w1_u, length.out = k)
  seq_w1_t1 = seq(bounds_t1$w1_l, bounds_t1$w1_u, length.out = k)
  
  #vote rate for veterans/non-veterans across possible values of w1
  grid_bw_t0 = lapply(seq_w1_t0, function(w1) 
    weighted.mean( (bounds_t0$w0 +w1*x), n*(1-x) ) ) %>%
    unlist()
  grid_bb_t0 = lapply(seq_w1_t0, function(w1) 
    weighted.mean( bounds_t0$w0 + bounds_t0$c1 + bounds_t0$d1 * x + w1*(x - 1), n*x)) %>% 
    unlist()
  
  
  grid_bw_t1 = lapply(seq_w1_t1, function(w1) 
    weighted.mean( (bounds_t1$w0 +w1*x), n*(1-x) ) ) %>%
    unlist()
  grid_bb_t1 = lapply(seq_w1_t1, function(w1) 
    weighted.mean( bounds_t1$w0 + bounds_t1$c1 + bounds_t1$d1 * x + w1*(x - 1), n*x)) %>% 
    unlist()
  
  #Get first differences for veterans/non-veterans
  diff_bb = outer(grid_bb_t1, grid_bb_t0, `-`)
  diff_bw = outer(grid_bw_t1, grid_bw_t0, `-`)
  
  #Get second differences
  diff_mat =  diff_bb - diff_bw
  
  #Get heterogeneous effects
  r_bb_t0 = lapply(seq_w1_t0, function(w1) 
    bounds_t0$w0 + bounds_t0$c1 + bounds_t0$d1 * range(x) + w1*(range(x) - 1)) %>% 
    #lapply(function(b) b[2] - b[1]) %>%
    do.call(., what = "rbind")
  
  r_bw_t0 = lapply(seq_w1_t0, function(w1) 
    (bounds_t0$w0 +w1*range(x))) %>%
    #lapply(function(b) b[2] - b[1]) %>%
    do.call(., what = "rbind")
  
  r_bb_t1 = lapply(seq_w1_t1, function(w1) 
    bounds_t1$w0 + bounds_t1$c1 + bounds_t1$d1 * range(x) + w1*(range(x) - 1)) %>% 
    #lapply(function(b) b[2] - b[1]) %>%
    do.call(., what = "rbind")
  
  r_bw_t1 = lapply(seq_w1_t1, function(w1) 
    (bounds_t1$w0 +w1*range(x))) %>%
    #lapply(function(b) b[2] - b[1]) %>%
    do.call(., what = "rbind")
  
  #Get differences in voting rates for high - low x
  diff_t0 = r_bb_t0[,2] - r_bb_t0[,1]
  diff_t1 = r_bb_t1[,2] - r_bb_t1[,1]
  
  diff_t0_w = r_bw_t0[,2] - r_bw_t0[,1]
  diff_t1_w = r_bw_t1[,2] - r_bw_t1[,1]
  
  #Differential time trend for high/low x.
  het_bb = outer(diff_t1, diff_t0, `-`)
  het_bw = outer(diff_t1_w, diff_t0_w, `-`)
  
  #
  #minimum/maximum values of
  #- AVERAGE DRIFT for veterans;non-veterans
  #- HETEROGENEOUS DRIFT for veterans;non-veterans
  #
  min_eff_v_l = outer(r_bb_t1[,1], r_bb_t0[,1], `-`)[diff_mat > 0] %>% min
  min_eff_v_u = outer(r_bb_t1[,2], r_bb_t0[,2], `-`)[diff_mat > 0] %>% min
  max_het_v = het_bb[diff_mat > 0] %>% max
  max_eff_nv_l = outer(r_bw_t1[,1], r_bw_t0[,1], `-`)[diff_mat > 0] %>% max
  max_eff_nv_u = outer(r_bw_t1[,2], r_bw_t0[,2], `-`)[diff_mat > 0] %>% max
  max_het_nv = het_bw[diff_mat > 0] %>% max
  
  #ASSUMING NO NEGATIVE EFFECTS
  #no_neg_effects = ((outer(r_bb_t1[,2], r_bb_t0[,2], `-`) > 0) & 
   #                   (outer(r_bb_t1[,1], r_bb_t0[,1], `-`) > 0)) %>% which(arr.ind = T)
  scaleFUN <- function(x) sprintf("%0.1f", x)
  plot_bb = data.table(v_diff_l = outer(r_bb_t1[,1], r_bb_t0[,1], `-`) %>% as.vector,
                       v_diff_u = outer(r_bb_t1[,2], r_bb_t0[,2], `-`) %>% as.vector,
                       ate = diff_mat %>% as.vector)
  v_plot = ggplot(plot_bb, aes(x = v_diff_l, y = v_diff_u, colour = ate)) + 
    geom_point() + theme_bw() + 
    scale_color_gradient2("ATT", low = 'red', high = 'blue',
                          midpoint = 0, breaks = seq(-.2,.4,.2)) +
    geom_hline(yintercept = plot_bb[abs(ate) == min(abs(ate)), v_diff_u]) +
    geom_vline(xintercept = plot_bb[abs(ate) == min(abs(ate)), v_diff_l]) +
    geom_hline(yintercept = plot_bb[abs(ate-att) == min(abs(ate-att)), v_diff_u], color = 'red') +
    geom_vline(xintercept = plot_bb[abs(ate-att) == min(abs(ate-att)), v_diff_l], color = 'red') +
    xlab("Mean Change in Pro-Suffrage Vote (low enlistment)") +
    ylab("Mean Change in Pro-Suffrage Vote (high enlistment)")
  
  plot_bw = data.table(v_diff_l = outer(r_bw_t1[,1], r_bw_t0[,1], `-`) %>% as.vector,
                       v_diff_u = outer(r_bw_t1[,2], r_bw_t0[,2], `-`) %>% as.vector,
                       ate = diff_mat %>% as.vector)
  nv_plot = ggplot(plot_bw, aes(x = v_diff_l, y = v_diff_u, colour = ate)) + 
    geom_point() + theme_bw() + 
    scale_color_gradient2("ATT", low = 'red', high = 'blue',
                          midpoint = 0, breaks = seq(-.2,.4,.2)) +
    geom_hline(yintercept = plot_bw[abs(ate) == min(abs(ate)), v_diff_u]) +
    geom_vline(xintercept = plot_bw[abs(ate) == min(abs(ate)), v_diff_l]) +
    xlab("Mean Change in Pro-Suffrage Vote (low enlistment)") +
    ylab("Mean Change in Pro-Suffrage Vote (high enlistment)") + 
    geom_hline(yintercept = plot_bw[abs(ate-att) == min(abs(ate-att)), v_diff_u], color = 'red') +
    geom_vline(xintercept = plot_bw[abs(ate-att) == min(abs(ate-att)), v_diff_l], color = 'red') +
    scale_y_continuous(labels = scaleFUN)
  
  out = list(v_plot = v_plot, nv_plot = nv_plot, 
         min_eff_v_l = min_eff_v_l ,
         min_eff_v_u = min_eff_v_u,
         max_het_v = max_het_v ,
         max_eff_nv_l = max_eff_nv_l,
         max_eff_nv_u = max_eff_nv_u ,
         max_het_nv = max_het_nv 
  ) 
  return(out)
}